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Abstract 

We introduce a new method of obtaining guaranteed enclosures of the 
eigenvalues of a variety of self-adjoint differential and difference operators 
with discrete spectrum. The method is based upon subdividing the region 
into a number of simpler regions for which eigenvalue enclosures are already 
available. 
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1 Introduction 

A rigorous method of obtaining enclosures of the eigenvalues of self-adjoint oper- 
ators has recently been described by Goerisch and Plum ||, |7|, |], |§. It depends 
upon having a soluble comparison operator, from which a controlled homotopy is 
carried out. In this paper we introduce a new method which has the advantage of 
not requiring such a comparison operator, and apply it to a variety of examples. 
In Sections 2 to 5 we consider Sturm-Liouville operators in some detail. Sections 6 
and 7 describe how to adapt the method to higher order operators and systems, still 
in one dimension. In Sections 8 and 9 we treat discrete Laplacians on graphs, while 
in Section 10 we consider the Laplacian acting in a bounded region in Euclidean 



space. The method can be applied to second and higher order elliptic differential 
operators with variable coefficients, but we do not present the details here. 

We distinguish between computing an eigenvalue in floating point arithmetic, and 
obtaining guaranteed enclosures. When using the word 'enclosure' we shall always 
understand that the calculation is mathematically rigorous, and that the compu- 
tations are done in interval arithmetic. Most numerical computations do not give 
proofs that the values obtained are correct, but depend upon the experience of the 
person who writes or uses the program concerning its range of reliability. With 
guaranteed enclosures on the other hand, the value obtained is known to be correct 
within the stated error bounds, unless there is an actual error at some stage of the 
computation. 

There are already several methods of computing the eigenvalues of a Sturm-Liouville 
operator H acting in L 2 (a,/3), and higher order analogues. The most obvious 
one, called shooting, solves the initial value problem for the eigenvalue equation 
Hf = Xf subject to the given boundary conditions at a and then varies A until 
the boundary condition at (3 is also valid. Most shooting programs do not try to 
give guaranteed error bounds. Although this is entirely possible ||, the method is 
difficult to implement numerically if the potential is singular at both ends of the 
interval. 

A second method, introduced by Goerisch || and Plum 0, §, §, obtains guar- 
anteed enclosures on the eigenvalues of a self-adjoint operator H by a continuous 
homotopy method, starting from a simpler operator. This is often exactly soluble, 
but a minimum requirement is that one can obtain sufficiently good rigorous lower 
bounds on its eigenvalues. Our method is similar to theirs in that it also uses a 
homotopy from a simpler operator. However they consider a continuous homotopy 
in some parameter, which often changes the coefficients smoothly to those of an 
operator with constant coefficients, while we consider a discrete homotopy in cer- 
tain internal boundary conditions which we choose to insert. We have compared 
our variation of the homotopy method with theirs for some of the examples Plum 
solves, and it appears to be substantially more efficient. In higher dimensions we 
are able to treat examples which are beyond the earlier method, because of the 
non-existence of an exactly soluble operator possessing a continuous homotopy to 
the given operator. 

One may obtain rigorous upper bounds on any specified number of eigenvalues by 
means of the Rayleigh-Ritz (RR) or variational method ]IJ. The starting point 
is the determination of accurate approximations to the eigenfunctions by a non- 
rigorous auxiliary calculation, possibly an inverse iteration method. Once these 
have been found one starts again using RR to obtain rigorous upper bounds on the 
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eigenvalues of the operator H in interval arithmetic. 



The lower bound is obtained by the method of Temple- Lehmann (TL) which also 
depends upon the choice of suitable test functions 0, |9|, [T0| . However, in this 



case one also needs to have crude lower bounds on the eigenvalues, and these are 
precisely what is missing at the rigorous level. More precisely if the eigenvalues of 
H, written in increasing order and repeated according to multiplicity, are {A n }^L , 
then in order to obtain an accurate lower bound on A n for some n using TL one 
needs already to be in possession of a number p such that 

A n < P < A n+ i 

where p is not too close to A n . There are three possible methods of obtaining such 
lower bounds. 

(i) One might hope that the upper bound on A n+ i is fairly accurate and take p to 
be a slightly smaller number. This idea cannot be turned into a rigorous procedure 
and will not be discussed further. 

(ii) One can use the Goerisch-Plum coefficient homotopy method of obtaining 
enclosures for many operators. 

(iii) One can use a boundary condition homotopy method. The description of this 
new method is the main contribution of this paper. 

In the above we have not mentioned the extra complications which arise if A n is 
degenerate or nearly so. There are well-known modifications of TL which deal with 
this problem fl], [7|, PH|, but we did not want to over-complicate the discussion at 
this stage. 



2 Neumann decoupling 

Let if be a Sturm-Liouville operator acting in L 2 (a,j3). We assume that H is of 
the form 

Hf{x) '■ = ~-L{ a{x) %} +v{x)nx) 

where a is a positive function in C l [a,f3] and V G L°°[a, /?]. We assume Neumann 
boundary conditions (NBC) in order to emphasise that the method does not depend 
upon the very strong monotonicity properties which hold for Dirichlet boundary 
conditions 

Our method is based upon decoupling the interval (a, 0) into 2 N subintervals; in 
most examples considered by the author one can take N = 3 or N = A. The 
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subintervals do not need to be of equal length, but this is the easiest choice to 
make. We put 

a = a® < a.\ < . . . < a 2 N = (3 

and let Hi denote the restriction of H to L 2 (aj_i, ctj) subject to NBC. We then 
define An to be the sum of the Hi, so that An once again acts in L 2 (a,P). The 
operators H and An have the same quadratic form 

Q(f) := [\a(x)\f'(x)\ 2 + V(x)\f(x)\ 2 }dx 

but with different quadratic form domains Q{H) and Q(An)- The space of all C 1 
functions on [a, (3] is a quadratic form core for H, but to obtain a quadratic form 
core for An one must allow the functions to have arbitrary jump discontinuities at 
each aj. Since Q(H) C Q(A N ), the RR method [|IJ shows that the eigenvalues of 
A^ are less than or equal to those of H . 

We define intermediate operators A n acting in L 2 (a, P) for < n < N, by a similar 
method. The component operators of A n are similar to those of An but only using 
the points ojj where i = j.2 N ~ r for < j < 2 r . At each stage the quadratic form 
domain decreases, leading to the operator inequalities 

A N < A N -i < ■ ■ ■ < A) = H. 

The following theorem enables the eigenvalues of A n to be computed rigorously 
using TL once one knows those of A n+ i. Since the passage from A n+ \ to A n 
consists of putting together intervals in pairs, and the set of eigenvalues of A n is 
simply the collection of all eigenvalues of its component operators Hi, it is sufficient 
to deal with the following special case, which also describes the passage from A\ 
to A = H. 

Let a < 7 < (3 and let {A^}, {fii}, {^t}, denote the eigenvalues of the operators 
Hi, H 2 and H associated with the intervals (0,7), (7,/?), (a, P) respectively, all 
subject to NBC. Finally let 

:= {\i} U {^} 

subject to re-ordering in increasing order and repeating according to multiplicities. 
Then {o"j} are the eigenvalues of the operator A := H\ + H 2 . 

Theorem 1 The eigenvalues {a{\ and {z/j} interlace in the sense that 

Oi<Vi< a i+ i 

for all i. Moreover these are strict inequalities unless the derivative of the relevant 
eigenfunction of H vanishes at the point 7. 
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Proof The idea is that H differs from A by a rank one perturbation in a certain 
singular sense. More precisely let s < <Jq and compare (A + s) -1 with (H + s) _1 . 
Both have Green functions which can be computed from the two fundamental 
solutions of the differential equation 



If one computes the difference of the two kernels one finds that it is a rank one 
operator. The inequality which we want is equivalent to 



and this holds whenever one has a positive rank one perturbation, by an application 
of the min-max principle. For an alternative proof see Theorem 5. 

The eigenvalues of H\ are all distinct, as are the eigenvalues of H 2 , because Hi are 
Sturm-Liouville operators. However there may be coincidences between the two 
sets of eigenvalues, which imply that <7j = (T i+1 . This is not a serious problem, 
but it can usually be avoided if desired by moving the point 7 slightly. Assuming 
that this has been done it is not possible that Oi — i/f. this equality would imply 
that the corresponding eigenfunction of H happens to have zero derivative at 7, in 
which case it is also an eigenfunction for both H\ and H 2 , so <Tj = 

3 The Enclosure Algorithm 

We start with a subdivision of (a, f3) into 2 N parts which is fine enough for us 
to be able to obtain disjoint enclosures on the eigenvalues of each component Hi 
by comparison with constant coefficient operators. From here onwards we suppose 
that we are only interested in obtaining enclosures on those eigenvalues of H which 
are less than some pre-assigned number E. The larger the value of E, the larger one 
must take N in order to be able to start the procedure. The following lemma shows 
that it is sufficient for the coefficients to be close to constant in each interval. We 
only consider the case of H itself for notational simplicity, but the lemma should 
actually be applied to each component Hi of A^. 

Lemma 2 Suppose that ao < a(x) < a\ and vq < V(x) < v\ for all x e (cx,f3). 
Then the eigenvalues {\i} of H satisfy 




(<Ti + S)' 1 > {Vi + S)' 1 > (<7 i+ i + S) 



-1 



a ir 2 i 2 /(/3 - a) 2 + v < A* < a^i 2 /^ - a) 2 + v ± . 
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These enclosure intervals are disjoint for all eigenvalues less than a given number 
E' if 

< V! - v < a 7r 2 /{(3 - a) 2 

and 

0<v 1 -v < M 2 a n 2 /(f3 - a) 2 — (M — l) 2 ai ir 2 /(P - a) 2 
where M is the smallest integer such that 

E' < v + M 2 a 7r 2 /([3 - a) 2 . 

Proof The first inequality follows by comparing H with the obvious constant co- 
efficient operators, whose eigenvalues are exactly computable. The proof of the 
second uses the observation that 

v o + m 2 a vr 2 / {(5 — a) 2 — v x — {m — l) 2 ai-n 2 / (j3 — a) 2 

is a concave function of m which must therefore take its minimum value on the 
interval [1, M] at one of its ends. 

Note Since the initial intervals (aj_i,a:j) are quite short, the integer M in the 
above lemma may be quite small and the above lemma may not impose strong 
conditions on the constants v , v±, a , a\. 

The algorithm for obtaining enclosures of the eigenvalues of H has several stages: 

Stage 1 We have to choose an initial subdivision of the interval such that 

each of the subintervals (ctj_i, «j) satisfies the conditions of Lemma 2. This can be 
done in several ways and is discussed further in Section 4. 

Stage 2 We choose a number E' > E, for example E' := 9E/8, and put E n : = 
E + n(E' - E)/N for all < n < N + 1. 

Stage 3 We subdivide (a, (3) as described above for a value of iV which is large 
enough for us to obtain disjoint intervals which enclose each of the eigenvalues of 
each component Hi of A N up the number E N+ i. 

Stage 4 We use RRTL to obtain accurate enclosures of each of the eigenvalues 
of each component Hi up to the number En- Putting these together in pairs we 
obtain rough enclosures of the eigenvalues of each component Hj of An-\ by virtue 
of Theorem 1. These enclosure intervals overlap very slightly because the previous 
accurate enclosures were not perfect. 
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Stage 5 We apply RR to each component Hj of A N -i to obtain smaller upper 
bounds on each of the eigenvalues of each Hj and so to convert the above into 
rough but nevertheless disjoint enclosures of the eigenvalues of each Hj up to En- 

Stage 6 We apply TL to obtain accurate enclosures of each of the eigenvalues of 
each component Hj of A N -i up to E N -i- 

Stage 7 We repeat the process inductively until we reach accurate enclosures of 
each of the eigenvalues of H up to E. 

Some comments are in order. 

The introduction of the sequence E n at Stage 1 is needed because TL requires a 
significant gap above any eigenvalue to be estimated. If there is an eigenvalue very 
close to the upper limit E n at any stage then that eigenvalue cannot be estimated 
accurately. 

When the eigenvalues af two adjacent operators are combined in Stage 3 it may 
happen that two eigenvalues of the new list created coincide to a high degree of 
accuracy. This is one possible cause of the problem mention in the next paragraph. 

The procedure in Stage 5 may occasionally fail because RR may not decrease 
the upper bound on an eigenvalue enough to make the intervals disjoint. This is 
handled by using a higher order version of TL whenever this occurs. In principle 
this could occur for all eigenvalues, in which case the algorithm might halt, but 
this is extremely unlikely unless there is a symmetry of the underlying problem, 
which should have been taken into account before starting the computation. 

Ultimately we do not guarantee either that the algorithm finishes or that the 
results which it yields are of the desired accuracy, but only that if the algorithm 
does finish then the enclosures obtained are correct. If the enclosures are not 
sufficiently accurate, then one must start again with a larger test function space. 

Although we specified that the second order coefficients a(x) of the differential 
operator should be C 1 , there is no difficulty in accommodating simple jump dis- 
continuities. Once one has determined the location of these points, they should be 
included in the partition {ai}f =0 of the interval (a, f3). The discontinuity of a(x) 
at a point 7 imposes an effective internal boundary condition on H at 7, which 
must be taken into account when specifying its operator domain, but has no effect 
on its quadratic form domain. 

One way of estimating the total computational effort is to count the number of 
distinct operators for which we have to compute some of the eigenvalues accurately. 
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At the level n this is 2 n , so the total number is 2 N+1 — 1. 

Since parallel machines will become more important, it should be noted that the 
computations of the eigenvalues of the different operators Hj at any particular 
level are entirely independent, and may be carried out simultaneously. Thus on 
a parallel machine the total computational effort is proportional to N + 1. In all 
the examples which we have considered this means the algorithm has only three 
or four steps! 

Both of the above estimates of computational effort are too pessimistic. At the 
higher levels it may be seen in the examples we analyse below that the number of 
eigenvalues of each operator to be computed is very small, because the eigenvalues 
are far apart. So the computation is much faster at the higher levels than indicated 
above, whether or not one has a parallel machine. 

It is clear that the same procedure may be used irrespective of the boundary 
conditions at a and /3. It may also be applied to potentials which are singular at 
the end points provided one has crude bounds on the eigenvalues to replace those 
of Lemma 2. Its extension to systems and to higher order differential operators is 
described in Sections 6 and 7. 

4 The Subdivision of (a,/3) 

We have suggested above that the subdivision of (ee,/3) should be defined by 

Oii := a + i(P - a) /2 N 

for < % < 2 N , where the size of N is determined by Lemma 2 as indicated in 
the algorithm. However it is possible that when combining two eigenvalue lists 
one finds that two eigenvalues coincide or are undesirably close. This is not an 
insuperable problem since one can use a higher order version of TL to obtain the 
required lower bounds on the eigenvalues. However, there is a systematic way of 
avoiding it, unless one of the eigenfunctions has an interval of constancy. 

We first emphasise that there is no hope of obtaining accurate enclosures of the 
eigenvalues of H unless there is some other non-rigorous method of computing 
the eigenvalues, such as unsupplemented RR or shooting, which in fact give good 
approximations to the eigenvalues and eigenfunctions. We use these computed 
eigenfunctions to choose the bisection point 7 of as described below. We 

then do the same for both of the subintervals (a, 7) and (7,/?) and so on until 
we have produced a fine enough subdivision of (a,/3) according to the criterion 
of Lemma 2. If we have misled ourselves about the best choice of the points ctj 
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then nothing is lost, because we can still use the above algorithm. If however, the 
approximations to the eigenfunctions are accurate enough then the method we now 
describe will have prevented the problem mentioned above. 

Let us suppose that there are k + 1 eigenvalues of H less than E, and that the 
corresponding eigenfunctions f r have zero derivatives at p(r) points for each < 
r < k. Then we choose 7 somewhere near the centre of (a,/?) but not at or near 
to any of the above points. Since there are P := p(0) + . . . + p(k) such points 
altogether there exists 7 G (3a/4 + /?/4, a/4 + 3/9/4) which is at a distance at least 
(f3 — a)/AP from each of the points. 

Whether or not it is worth using this iterative procedure for selecting the subdivi- 
sion of (at, [3) remains to be seen. In the two cases solved below, it appears that 
using a uniform subdivision is perfectly satisfactory. 

There is an entirely different reason for choosing a non-uniform subdivision of 
(a,P), if the coefficients of H vary substantial from one part of the interval to 
another. If the potential is bigger than the number E in some interval, then one 
should make that entire interval one of the (aj_i,Oj), however big it is, because 
there will be no relevant eigenvalues associated with it. More generally the size of 
each interval (a^i,^) should be as big as possible subject to being able to ob- 
tain disjoint enclosures of all of the eigenvalues of Hi less than E. This procedure 
reduces the number of operators for which one has to compute some of the eigen- 
values. A more thorough investigation might involve the uncertainty principle, but 
Lemma 2 suffices for most purposes. 

5 Examples 

We illustrate our general theory with two numerical examples, which are solved 
using shooting and floating point arithmetic, not using interval arithmetic as is 
actually required. There are two reasons for this, the first being that our goal here is 
only to examine the feasibility of the method, not to create a new software package. 
The second is that one should not use a high-powered technique for obtaining 
eigenvalue enclosures until one has a good idea of the approximate location of the 
eigenvalues. This information cannot be used in the final computation because it 
is not rigorous, but it may indicate problems which need special attention in the 
rigorous computation. 

The two examples were studied in detail by Plum |7j] using a continuous homotopy 
in the coefficients. 
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Example 3 Let H be the operator defined by 

j2 £ 

Hf(x):=-^L + 8cos(x) 2 f(x) 

acting in L 2 (0,7r) subject to NBC. Plum obtained enclosures on the eigenvalues 
ranging from 

fi = 2.48604311^ 

to 

/i 8 = 68.03175^ 

using his homotopy method, RRTL and interval arithmetic. Putting N = 2 the 
conditions of Lemma 2 are satisfied for any choice of E for each of the components 
Hi, 1 < i < 4, with vi — vq = 4, ao = a± = 1 and oti — a^i = it/4. Now let K±, K 2 
be the two operators at level one, acting in the intervals (0, it/ 2) and (tt/2, it). We 
have computed the eigenvalues of all of the operators above up to the limit E — 70. 

For % — 1 , 4 Lemma 2 yields the crude enclosures 

4 < < 8, 20 < fa < 24, 68 < fi 2 < 72 

while more accurate, but non-rigorous, calculations provide 

Ho ~ 6.454, /ix ~ 22.450, ^ 2 ~ 70.515 

For ? = 2,3 Lemma 2 yields the crude enclosures 

< no < 4, 16 < nx < 20, 64 < /i 2 < 68, 144 < /i 3 

while more accurate calculations yield 

Ho ~ 1.364, /ix ~ 17.693, /i 2 ~ 65.503 

We now join together the eigenvalue lists of Hi and H 2 to obtain the list 

1.364, 6.454, 17.693, 22.450, 65.503, 70.515 

If these values are indeed accurate then according to Theorem 1 they interlace 
the eigenvalues of Ki and provide the basis for the use of RRTL for obtaining 
accurate enclosures of the eigenvalues of K\. The eigenvalues of Ki are (again 
non-rigorously) 

2.486, 9.173, 20.141, 40.057, 68.032 

These coincide with the eigenvalues of K 2 since we have not made use of the 
symmetry of the operator about x = tt/2. When we combine this list with a 
second copy of itself the resulting list interlaces the eigenvalues of H, namely 

2.486, 6.397, 9.173, 13.370, 20.141, 29.084, 40.057, 53.042, 68.032 
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It would have been possible to avoid the coincidence of the eigenvalues of K\ and 
K 2 by starting with the partition 0, 0.6, 1.2, 2.1, tt instead of the partition into 
equal subintervals. 

The above computation involves determining the eigenvalues of 7 operators at 3 
different levels. At the top level we only had to compute the first 3 eigenvalues of 
each operator Hi, at the middle level we had to compute 5 and at the bottom level 
we had to compute 9. 

By comparison Plum computed the eigenvalues of 10 intermediate operators, with 
less scope for parallelization since each computation depended on the previous one. 
We computed a total of 31 eigenvalues, while Plum computed at least 90. Plum 
was not, however, particularly concerned with minimising numerical effort in his 
paper. 

Example 4 We consider the operator 

d 2 f 

Hf(x) := — r4 + I000xf(x) 
dx z 

acting on L 2 (0, 1) subject to DBC This is essentially the same as Example 2 of 
Plum H|, who obtained enclosures on the eigenvalues ranging from 

Ho = 233.8107^ 

to 

/i 9 = 1508.10?! 

The unpublished enclosures of Lohner || , obtained by shooting, are considerably 
more accurate. We put N := 3, «j := i/8 for < i < 8 and E := 1000. A 
modification of Lemma 2 to cope with the DBC at 0, 1 yields the following crude 
eigenvalue enclosures. 

In (0, 1/8) we have initially 

157 < fio < 283, 1421 < Hi < 1547 
and then more accurately 

Ho ~ 245.225, Hi - 1486.798 
In (1/8, 1/4) we have initially 

125 < Ho < 250, 756 < Hi < 882 
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and then more accurately 



fio ~ 185.471, nx ~ 820.761 

We omit the results for the other intervals, each of which involves the computation 
of only two eigenvalues, all higher eigenvalues being bigger than 1500. We now 
consider level 2. Putting the previous lists together in pairs and considering the 
interval (0, 1/4) we obtain the crude bounds 

185.471 < < 245.225 < \i x < 820.761 < fi 2 < 1486.798 < fi 3 

In fact the upper bound on each /ij is slightly bigger than the lower bound on 
jtii+i, because the number separating them is not exact, but the use of RR reduces 
the upper bound on each eigenvalue substantially and so yields disjoint enclosures 
of the eigenvalues (or actually would do so if the calculations were rigorous). All 
higher eigenvalues are greater than 1500. The accurate eigenvalues for the interval 
(0,1/4) are 

li = 205.942, //! = 490.938, /i 2 = 1115.419 

We omit the further computations at levels 2 and 1. The full list of eigenvalues of 
Ai is 

233.705, 400.348, 532.152, 601.881, 748.110, 825.999, 1007.897 
and interlaces the list of eigenvalues of H, namely: 

233.811, 408.795, 552.056, 678.679, 794.738, 906.461 

which agree with the enclosures of Plum [0. We next comment on the amount of 
computation needed by our method. 

The total number of operators considered by our method, is 15, compared with 50 
in Plum's method, since he puts 5 = 0.02. If one has a parallel machine then the 
relevant quantity is the number of levels, namely 4. For each operator at level 3 we 
needed to compute 1 or 2 eigenvalues. For each operator at level 2 we computed 
3 eigenvalues. For the two operators at level 1 we computed 5 and 3 eigenvalues. 
Finally we computed all 6 eigenvalues of H in the interval [0, 1000], making a total 
of at most 42 eigenvalues computed. Plum's method involves the computation of 
at least 300 eigenvalues, but he did not attempt to minimise this number. 

6 Higher Order Operators 

The procedure which we described above can be modified to treat higher order 
differential operators in one dimension. The difference in the higher order case is 
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that decoupling an interval into two parts by introducing a Neumann boundary 
condition is not equivalent to a rank one perturbation. However it is still of finite 
rank, as we will now explain. 



Let H be defined formally on L 2 



(a, (3) by 



Hf(x) := (-1) 



rn 



dx m 



a(x) 



dx m 



d m f 



I 



where a G C m [a,j3] is positive. Our method can also deal with more complicated 
operators involving lower order terms. We assume Neumann boundary conditions, 
in the sense that we take the quadratic form of the operator to be 



with domain the Sobolev space W m ' 2 (a, (3). It is known that Q is closed on this 
domain, and we define H to be the non-negative self-adjoint operator associated 
with the form in the standard manner. 

Functions in W m ' 2 (a, (3) are continuous on [a, j3] along with all derivatives of order 
less than m. Given a < 7 < (3 we introduce a Neumann boundary condition at 
7 by replacing W m,2 (a, (3) by the space Q m in which we allow the functions and 
their first m — 1 derivatives to have simple jump discontinuities at 7. Let H m be 
the corresponding operator on L 2 (a, f3). We now define a chain of operators H r for 
< r < m with H := H . Each of them is associated with the same form Q but 
on different domains Q r . We define Q r to be the space of functions in Q m such 
that all derivatives of / from the order r to m — 1 inclusive are continuous at 7. 
Thus Q r C Qr+\ for all r, each being of co-dimension one in the next. 

Theorem 5 Let H and K be two non-negative self-adjoint operators on a Hilbert 
space 7i such that their quadratic forms coincide on their common domain. Suppose 
also that the form domain Q{K) of K is a subspace of co-dimension 1 in Q(H). 
Finally suppose that H and K both have purely discrete spectrum and that their 
eigenvalues written in increasing order and repeated according to multiplicity are 
respectively {A n }^L and {/i„}^L . Then the two sets of eigenvalues interlace in 
the sense that 



Proof This is an immediate consequence of the min-max principle, since every 
subspace of dimension n of Q(H) is either already contained in Q(K) or intersects 
Q(K) in a subspace of dimension n — 1. 




dx 



for all n. 
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The application of this theorem to higher order operators is immediate. In order to 
remove a Neumann boundary condition at the point 7 we have to pass through a 
chain of operators H r with r decreasing from mtoO. At each stage the eigenvalues 
interlace, and this is the condition needed to apply the TL technique as described 
in Section 3. 

We have described the operator H r in terms of its quadratic form domain. This 
is sufficient for the application of the RR technique. However, the TL method 
depends upon the selection of test functions from the operator domain, so we 
need to describe this. We first comment that functions in any of the operator 
domains lie in C 2m ~ 1 [«, 7] + C 2 " 1-1 ^, f3] because of our smoothness assumption on 
the coefficient a(x). Also the weak derivative /( 2m ) (x) in each subinterval must lie in 
L 2 . Of course eigenf unctions are more regular and must lie in C 2m [a, 7] +C 2m [ 7 , (3\. 

We now specify the boundary conditions. The choice of quadratic form domain 
implies that if / lies in the operator domain then f^ r \x) = for x = a, (3 and for 
all m < r < 2m — 1, i.e. Neumann boundary conditions. We need to impose 2m 
boundary conditions at 7± to obtain a self-adjoint operator, and these are different 
for each operator H r . Our quadratic form assumption is that f^i'J—) = /^(7+) 
for all s such that r < s < m — 1. This corresponds to the assumption that 

/W( 7+ ) = /W( 7 -) for all r < s < m - 1 

( a /(m))(s)( 7+ ) = ( a /M)(*)( 7 _) f or a n < s < m - r - 1 
( a jF("»))W( 7 ±) = forallm-r<s<m-l 

for all / in the operator domain of H r , as one may see by carrying out some 
integrations by parts and requiring the boundary terms to vanish. 

The test functions chosen for the TL procedure must satisfy all of the above bound- 
ary conditions. One could use a space consisting of different polynomials in each 
subinterval, with the coefficents restricted to satisfy the boundary conditions at 
a, P, 7, but many other choices are possible. 

7 Systems of Ordinary Differential Equations 

A self-adjoint system of Sturm-Liouville operators is defined as an Operator H 
acting in L 2 ((a,(3), C m ) according to the formula 

We assume that a it j G C l [a, j3] and V it j G L°°[a, (5] for all i, j. We assume that both 
matrices are real symmetric for all x G [a, (3] and that a it j is uniformly positive 
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definite on [a, 0\. We finally assume that the operator satisfies NBC in the obvious 
sense for systems. 

The computation of the eigenvalues of H proceeds as in the scalar case with one 
exception. Namely the removal of an internal NBC involves a perturbation of 
rank m rather than of rank 1 as in the scalar case. We deal with this as we did 
for higher order Sturm-Liouville operators in the last section. Instead of writing 
out the details in the general case, we solve a simple example, which exhibits the 
essential features of the general case. 

Example 6 Put m = 2 and a>i,j{x) := Sij for all Let a < < (3 and let u, v 
be arbitrary non-negative numbers. Then define the matrix-valued potential V by 



V(x) 



u 



v v 

V V 



if a < x < 
if < x < (3. 



Let Hi, H 2 , H be the operators associated with the above expression acting in the 
intervals (ck,0), (0,f3), (a, (3) respectively, all subject to NBC. Let K be the 'same' 
operator acting in the interval (a, (3), subject to NBC at a, (3 and the following 
boundary conditions at 0, expressed in terms of the operator domain: 

/i(0+) = A(O-) 

/i(o+) = mo-) 

/a(0+) = 
/a(O-) = 0. 



If Ai is the operator Hi + H 2 acting in L 2 (a,[3), then the quadratic forms of 
Ai, K, H are all given by the expression 

W) := f{\fi\ 2 + \f'2\ 2 + E Vij(x)Mx)f&)}dx. 

Ai has the largest quadratic form domain, W 1 > 2 ((a, 0), C 2 ) + W^ 1 - 2 ((0, (3), C 2 ), while 
H has the smallest quadratic form domain W 1,2 ((a, (3), C 2 ), of codimension 2 in 
the previous one. In between these lies the quadratic form domain of K, which is 
the set of / G W l ' 2 ((a, 0), C 2 ) + W^' 2 ((0, (3), C 2 ) such that /i(0+) = /i(O-). 

Since each quadratic form domain is a subspace of codimension 1 of the previous 
one, the eigenvalues of the operators interlace in the sense of Theorem 5. The 
eigenvalues of Hi, H 2 are exactly computable, so these observations allow us to 
obtain enclosures of the eigenvalues of H using RRTL in the standard manner. 
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We have chosen this example because the eigenvalues of all four operators involved 
are essentially exactly computable, and it is easy to confirm the interlacing property 
directly. We put a = —1, (3 := 2, u = 2v = 100, and compute all of the eigenvalues 
of each operator up to E := 50. 

The eigenvalues of Hi consist of all numbers of the form u + n 2 Ti 2 /a 2 or m 2 iT 2 /a 2 , 
where m, n are non- negative integers. This yields the list: 

0, 9.870, 39.478, 88.826 

The eigenvalues of H 2 consist of all numbers of the form 2v + n 2 ir 2 /(3 2 or m 2 Tc 2 /(3 2 , 
where m, n are non- negative integers. This yields the list: 

0, 2.467, 9.870, 22.207, 39.478, 61.685, 88.826 

The eigenvalues of A\ are obtained by combining these two lists to obtain: 

0, 0, 2.467, 9.870, 9.870, 22.207, 39.478, 39.478, 61.685, 88.826 

The eigenfunctions of K and H are linear combinations of trigonometric and expo- 
nential functions, and the eigenvalues are obtained by solving certain trancendental 
equations associated with the boundary conditions at 0. The eigenvalues of K are 
approximately: 

0, 0.468, 4.298, 9.870, 12.288, 24.757, 39.478, 41.865, 63.639 
which interlace those of Ai. The eigenvalues of H are approximately: 

0.449, 1.609, 4.735, 11.746, 17.747, 27.360, 41.177 
which interlace those of K. 

It may be seen that although the eigenvalues do interlace as the theory predicts, 
the smallest eigenvalue of H is rather close to the second eigenvalue of K, a fact 
which does not help the efficiency of the TL method. The reason for this is that 
the coefficients u, v are rather large, and this has the effect of partially decoupling 
the two intervals. Accurate lower bounds on the smallest eigenvalue of H can be 
obtained by using a higher order version of the TL method. 

8 Operators on graphs 

The method which we have developed for Sturm-Liouville operators may be applied 
with modifications to elliptic partial differential operators and to discrete Lapla- 
cians on graphs. The first application demands the use of the quite complicated 
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machinery associated with the finite element method. We decribe here the second 
application, which is of independent interest, and also involves the theory of rank 
1 perturbations in certain situations. 

We define a graph to be a finite set X together with a set £ of directed edges. 
We assume that if e := (x,y) G £ then e := (y,x) G £. We define the associated 
Laplacian to be the operator acting on l 2 (X) with matrix 

-1 H(x,y)e£ 
deg(x) if x = y 
otherwise 

where deg(x) := G X : (x, y) G £} is the degree of x. 
The quadratic form corresponding to this matrix is 

QU)--=\ E \f(*)-f(y)\ 2 

(x,y)ee 

which is a non-negative Dirichlet form, with all of the structural consequences of 
this fact. We follow standard practice in referring to the eigenvalues of the operator 
A defined above as eigenvalues of the graph X. 

Although the matrix A is finite the determination of its eigenvalues is not straight- 
forward if the graph is very large, and one actually has the same problems in 
obtaining guaranteeed enclosures as for infinite-dimensional problems. The first 
main problem is the non-existence of standard comparison problems. There are 
very few finite graphs for which one can compute the eigenvalues exactly, and there 
is no possibility of using a change of variables to map a graph to a standard soluble 
one as in the case of partial differential operators. 

We present two procedures for obtaining enclosures of eigenvalues of finite graphs. 
The first applies to the case in which the graph is obtained from one for which 
one already has enclosures of the eigenvalues by the removal of a small number of 
chosen vertices. We leave the reader to formulate the corresponding lemma relating 
to the additional of a small number of vertices. 

Lemma 7 Let Y be a subset of X obtained by the removal of a small number of 
vertices, and let Q be the set of (undirected) edges of X which join points ofY to 
points of X\Y. Then one may compute enclosures of the eigenvalues of Y from 
those of X in #(G) homotopy steps. 

Proof Let {ej}™ =1 be some enumeration of the edges in Q. Let Ai be the matrix 
acting in l 2 (X) which is obtained from that of A by the removal of the contribution 
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of the edges ei, . . . , e«. Then 

A > A 1 > . . . > A n 

and each matrix is a rank 1 perturbation of the next one in the chain. It follows by 
an argument similar to that of Theorems 1 and 5 that the eigenvalues of A^ interlace 
those of A i+ i for every i. This is the fundamental requirement for transferring 
accurate enclosures of the eigenvalues from A { to A i+1 by the RRTL method. The 
operator A n is the direct sum of the discrete Laplacians of Y and X\Y. Since we 
have assumed that X\Y is small, its eigenvalues may be computed independently 
by a direct procedure. Removing these eigenvalues leaves those of Y. 

Let X be a finite subset of Z N for some N. The edges of X are defined to be those 
pairs x, y G X such that 

N 

^ ' | fir Ur\ 1 ■ 
r=l 

In this situation we have the bound deg(x) < 2N for all x G X. The operator A 
may be considered to be the discrete Laplacian on X subject to Neumann boundary 
conditions, since is always an eigenvalue of A, the corresponding eigenfunction 
being constant. The multiplicity of the eigenvalue equals the number of connected 
components of the graph. 

Example 8 Let X C Z 2 be the set 

{(m, n) : 1 < m < k, 1 < n < k} 
and let Y be obtained by the removal of the set 

Z := {(1,1), (1,2), (2,1)}. 

Then an enumeration of the (undirected) edges of Q is e\ := ((1,2), (2,2)), e 2 : = 
((2,1), (2, 2)), e 3 := ((1, 2), (1, 3)), e 4 := ((2, 1), (3, 1)). We chose k := 7 and 
computed the smallest 6 eigenvalues of the operators Ai. The interlacing property 
is verified. The eigenvalues of A^ coincide with those of the discrete Laplacian of 
Y, together with the eigenvalues 0, 1, 3 of the Laplacian of Z. The numbers in the 
table below are actually k 2 times the eigenvalues, so that they may be compared 
with the eigenvalues of —A on the unit square subject to NBC. 







Hi 


A*2 


1^3 


A*4 


fJ>5 


A 





9.705 


9.705 


19.410 


36.898 


36.898 


A, 





9.515 


9.705 


19.142 


33.499 


36.898 


A 2 





9.361 


9.574 


18.367 


30.187 


34.782 


A 3 





5.868 


9.540 


13.119 


23.836 


34.571 


A, 








9.095 


11.471 


23.049 


32.525 
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The interlacing property states that the number A immediately above any eigen- 
value \x of Ai in the table is a lower bound for the next eigenvalue of A$. Let us 
suppose that the values in the first row of the table are close to accurate enclosures 
of the eigenvalues of A, and that all of the other entries in the table are rigorous 
upper bounds, which we expect to be accurate. Using the interlacing property 
we deduce that n^(Ai) = 36.898. The fact that the (upper bounds on the) other 
eigenvalues of A 1 are widely separated enables us to use TL to confirm that they 
have been found accurately. Interlacing establishes that fi 5 (A 2 ) > 33.499, and we 
confirm that the other eigenvalues of A 2 have been computed accurately as before. 
If interval arithmetic has been used we end up with accurate enclosures of all of 
the eigenvalues of A4 up to and including /i 4 . The final reason why this procedure 
works is that the entries in the column labelled fi^ and the row labelled A4 decrease 
rapidly enough for TL to be an efficient method. 

The above example is purely illustrative: the same procedure can be carried out 
for values of k which are large enough for the problem of obtaining eigenvalue 
enclosures to be non-trivial. If at some stage an eigenvalue does not decrease 
enough from one stage to the next for TL, then either we have to proceed to a 
lower eigenvalue, or we must use a higher order version of TL. 

The above method is not suitable for obtaining enclosures of the eigenvalues of 
a graph which is far from any graph for which eigenvalue enclosures are already 
known. As a typical example we mention the set of all (m, n) G Z 2 which satisfy 
all three inequalities m 2 + n 2 < 16d 2 , (m — 2d) 2 + n 2 > d 2 and (m + 2d) 2 + n 2 > d 2 
where d is a large positive number. 

In cases such as the above we combine the continuous homotopy procedure intro- 
duced by Goerisch and Plum with the hierarchical homotopy method we intro- 
duced for Sturm-Liouville operators. The idea is to subdivide X into several more 
or less convex parts each of which is small enough that eigenvalue enclosures can 
be obtained by a direct method. These parts are then joined together in pairs as 
described below, obtaining eigenvalue enclosures for the larger parts. If the initial 
subdivision is into k := 2 N parts, then after the first stage one has 2 N ~ 1 parts, 
and the procedure terminates after N stages. It remains to describe how to join 
together two subgraphs. 

Let X = Y U Z where Y, Z are disjoint subgraphs, let Q be the set of edges joining 
points of Y and Z, and let T be the complement of Q in the set £ of all edges of 
X. Given < s < 1, let A s be the matrix associated with the quadratic form 

Q.(f)-=l E l/(*)-/(y)l 2 + | E \f(x)-f(y)\ 2 - 

{x,y)er (x,y)eS 
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Lemma 9 The eigenvalues of A s are increasing real analytic functions of the pa- 
rameter s. The eigenvalue list of A is just the union of the two eigenvalue lists of 
Y and Z. At the other end A\ is the discrete Laplacian of X . 

Proof The first statement is part of received knowledge || , while the second de- 
pends upon the fact that A is the direct sum of the discrete Laplacians of Y and 
Z. 

The procedure for obtaining eigenvalue enclosures for X is similar to that of Go- 
erisch and Plum j|, [7|. We consider the operators A s ^ for a large enough chosen 
sequence = s < s± < . . . < s p — 1. If we have enclosures of the eigenvalues 
of A S (A then these provide lower bounds on the eigenvalues of A s n + x) which may 
be adequate to obtain enclosures of the eigenvalues of the latter operator by the 
RRTL procedure. Eigenvalue crossings may occur at certain values of s, but these 
are handled using the higher order TL procedure. 

Example 10 Let X := Y U Z C Z 2 where 

Y := {{m, n) : 1 < x < h — 1, 1 < y < x} 

Z := {(m,n) : h < x < 2h - 1, 1 < y < 2h - x} 

where h is some positive integer. Let the set £ of edges of X be those inherited 
from Z 2 as before. The undirected edges of Q are of the form (h — 1, r), (h, r) where 
1 < r < h — 1, and separate the triangle X into two smaller triangles. We list the 
7 smallest eigenvalues of A s below for h :— 8 and s := 0, 0.2, 1. A larger number 
of values of s were originally computed, but these are the only ones needed. 





/'o 


/'i 


M2 


V>3 


A*4 


A*5 












0.12061 


0.15224 


0.25330 


0.32139 


0.46791 


A .2 





0.04705 


0.13054 


0.23249 


0.27273 


0.42485 


0.49950 


A, 





0.07244 


0.13259 


0.27719 


0.33076 


0.51058 


0.60389 



Each eigenvalue fii is a monotonic increasing function of s. Suppose that we al- 
ready know that the numbers in the first row are accurate approximations to the 
eigenvalues of Aq, and that the other numbers are rigorous upper bounds to the 
corresponding eigenvalues, as determined by RR. We use the monotonicity to de- 
duce that p,e(Ao,2) > 0.46791. Using TL we then confirm that ^5(^0.2) is accurate. 
Monotonicity implies that ^(Ai) > 0.42485 and we are finally able to confirm the 
accuracy of Hj(Ai) for j = 4, 3, 2, 1, in turn by TL. If all of the computations have 
been done in interval arithmetic, we have obtained enclosures of the eigenvalues of 
A 1 from those of A . 
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The values of (J,i(A s ) for s = 0.. 1(0.1) are 

0, 0.03153, 0.04705, 0.05559, 0.06085, 0.06439 
0.06693, 0.06882, 0.07030, 0.07148, 0.07244 

This list exhibits a common pattern of rapid increase for small values of s followed 
by little change for large values. In fact it follows from RR that yt,\{A s ) is a concave 
function of s, but this need not be true for higher eigenvalues. 

A precondition for applying the above method is that eigenvalue enclosures of the 
parts Xi, . . . ,Xk into which we subdivide X should already be known. This may 
be achieved by making each Xi small enough so that all of its eigenvalues can be 
computed by a direct method. If some of the parts are rectangles or one of a very 
small number of other graphs then their eigenvalues may be exactly known. The 
following argument shows that for some purposes we may dispense with knowledge 
of accurate enclosures of the eigenvalues of the parts entirely. It requires instead 
an assumption about the geometry of the parts, expressed initially in terms of a 
lower bound on their first non-zero eigenvalues. 

Given a > we say that a finite graph (X, £) lies in C a if its first non-zero eigenvalue 
Hi satisfies 

fi! > ad(X,£)~ 2 

where d(X, £) is the diameter of the graph. We investigate the geometric signifi- 
cance of this condition in the next section. The value of b in the following theorem 
indicates how small the individual subsets in a partition of X need to be in order 
to be able to obtain accurate enclosures of the eigenvalues of X without already 
possessing accurate enclosures of the eigenvalues of the subsets. 

Theorem 11 Suppose that a > and that {Xj}^ =1 is a partition of the graph X , 
each subset of which lies in C a . Suppose also that b > and that 

d(Xi,£i) < U(X,£) 
b 

for all 1 < i < k. Then one may obtain accurate enclosures of all eigenvalues of 
X less than 

E := ab 2 d(X,£)- 2 
by a continuous homotopy method. 

Proof The conditions of the theorem imply that the first non-zero eigenvalue of 
each part is at least as big as E. Let (Y, JF) be the union of two of the subgraphs 
(X i: £i), and let (Y, T) be obtained by including also those edges of the graph 
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(X, £) which connect the two parts of Y . The spectrum of (X, J 7 ) below E consists 
of the eigenvalue with multiplicity 2 and nothing else. This precise if unusual 
information is enough to obtain accurate enclosures of the spectrum of (Y, 7) 
below E by the Goerisch-Plum continuous homotopy method. By a repetition of 
this method, carried out in a hierarchical manner, one eventually obtains accurate 
enclosures of the spectrum of (X, £) below E. 

9 The Poincare Inequality 

In order to implement the above ideas one needs to obtain geometric conditions 
which imply that the first non-zero eigenvalue of a connected graph (X, £) has a 
lower bound of order d~ 2 , where d is the diameter of the graph. Examples in || 
show that this is not always uniformly true for a family of graphs parametrised by 
the diameter as d — > oo, but their discrete version of the Poincare inequality can be 
rewritten to provide exactly what we need. We develop the theory of this section 
at a greater level of generality than before, because of its independent interest. 

Let b : £ — > (0, oo) be a positive weight function satisfying 6(e) = b(e) for all e G £. 
Let \X\ denote the number of points in X. Given a path 7 := (71, . . . ,7^), we 
define its length to be 

k 

l7l ^EKviJi)- 1 

i=l 

and then define the diameter d of X using this notion of length, in the usual 
manner. 

Let A be the operator on l 2 (X) associated with the quadratic form 

W):=~ E b(x,y)\f(x)-f(y)\ 2 . 
(x,y)ee 

It is easily seen that is an eigenvalue of A of multiplicity 1, the corresponding 
normalised eigenfunction satisfying 4>o(x) = |X| -1 / 2 for all x £ X. 

In order to obtain a lower bound on the first non-zero eigenvalue fii of A, we follow 
closely Diaconis and Stroock [|| (and Poincare). Suppose that T is a set of paths 
in X, one path from T joining every ordered pair of points x, y 6 X. We impose 
two constraints on the choice of this set of paths. The first, that 

\lx, y \ < ad 

for some a and all x, y G X, is self-explanatory. The second is that 

#{7 G T : e G 7} < Pd\X\ 
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for some (3 > and all e G £. Since the total number of paths in V is |X| 2 , this is 
a constraint on how well distributed the paths are. The assumption is in precisely 
the form needed for applications. 

Theorem 12 Under the above two assumptions we have 

1 



Proof If e := (x, y) G £ we put df(e) := /(y) — f(x). We have 

|x|||0-(0,0 o )0oir = \ E l0(*)-<Ky)l 2 

= ^ E I E ^(e)i 2 

< \ E l7x,»l E Ke) 150(e) | 2 

x,y€X ee-y x ,y 

< KQ(<P) 



where 



K : = sup E l7*,»l 

< ad#{j G T : e G 7} 

< a/3rf 2 |X|. 



The proof is completed by using the variational characterisation of n\\ 

W) 

||0-(0,0O)0O 112 



Upper bounds of a similar type on ^ are relatively easy to obtain by applying the 
variational inequality to suitable test functions, but we do not need them here. 

The following application of the above theorem is a typical building block for the 
implementation of the ideas in the last section. 

Theorem 13 Define X C Z 2 by 

X :={(i,j):l<i<n, l<j<f(i)} 
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where /{l, . . . , n} — > {1, . . . m} zs a non-decreasing function. Let A be the discrete 
Laplacian on X , corresponding to the choice 6=1 above. Then 

~~ (n + f (n) — 2) max{n, f (n)} ~ 
where d is the diameter of X . 



Proof It follows from the definition of X that d = n + f{n) — 2 and n < \X\ < 
nf(n). Our main task is to define the set T of paths. If i < %' then the path from 
(i,j) to (i',j') is the horizontal line from (i, j) to (i',j), followed by the vertical 
line from (i',j) to {i',j f ). Because / is monotonic, this is entirely contained in X. 
If i > i! we take a similar path. It may be seen that every path 7 has length at 
most n + f{n) — 2. The number of paths through any horizontal edge e G £ is 
at most n|X|, while the number through any vertical edge is at most f(n)\X\. A 
slight modification of the estimate of K in the last theorem completes the proof. 

The number of paths through any edge e G £ can be bounded more efficiently 
if further information about / is provided. If f(i) := 1 for 1 < i < n — 1 and 
f{n) := n, then the theorem yields /ii > n ^ 2 n-2) wm le one actually has /ii ~ ^ 
for large n. It would be valuable to determine the largest constant c such that 
fii > c/n 2 for all graphs of the type described in the above theorem, subject to 
f(n) < n. It appears that even the continuous analogue of this problem is unsolved. 



10 The Laplacian in N Dimensions 



We finally describe the modifications to the above ideas needed to provide enclo- 
sures of the eigenvalues of a partial differential operator. We will consider only 
the case in which H := A, acting in L 2 (Q) subject to NBC, where Q is a bounded 
region in R with piecewise smooth boundary. However, the same method applies 
to variable coefficient elliptic operators subject to other boundary conditions. By 
the eigenvalues of any region Q we mean the eigenvalues of —A acting in L 2 (Q) 
subject to NBC. 

If the region Q is diffeomorphic to the unit ball B Plum || has described a 
method of obtaining enclosures on the eigenvalues by transferring the operator 
to L 2 (B,m(x)dx) where m is a suitable positive weight, and then using his coef- 
ficient homotopy method. This cannot be adapted to treat the case in which O 
contains more than one hole. The method described below is capable of dealing 
with regions containing any number of holes. 
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Suppose we wish to find all of the eigenvalues of a region f2 which are smaller 
than a given number E > 0. The first step is to divide the region into subregions 
{Qi}i =1 each with a piecewise smooth boundary. We assume that enclosures of the 
eigenvalues less than E of each subregion are known, either because its eigenvalues 
are exactly computable or because it is small enough with a regular enough shape 
for to be its only eigenvalue below E. We also assume that the subregions can be 
recombined in pairs in a hierarchical fashion to recover the original set Q. The task 
therefore is to obtain enclosures of the eigenvalues of the union of two regions which 
have some common boundary when we already have enclosures of the eigenvalues 
of the individual regions. 

Let U, V be disjoint bounded connected regions in with piecewise smooth 
boundaries and suppose that their common boundary B is a non-empty (N — 1)- 
dimensional surface. Put Q := U U V U B . Suppose also that we have accurate 
enclosures of all the eigenvalues of each region up to the number E. If we combine 
the two lists of eigenvalues into a single increasing list {/ij(0)}, then this list pro- 
vides the eigenvalues of the operator H Q = —A acting in L 2 (fl) subject to NBC on 
dU U dV. Note that the number is an eigenvalue of multiplicity 2. 

We introduce a family of quadratic forms Q s defined for < s < oo, all having the 
same quadratic form domain T>q := W 1,2 (U) + W 1,2 (V). A core for this subspace 
consists of all functions on Q which are C 1 except that they are allowed to be 
discontinuous as one crosses B. Every function / 6 Z> nas L 2 boundary values on 
B which we denote by f± depending upon which side of B one approaches it from. 
We then define 

Qs(f) ■= I |v/| 2 + s / I/+-/-I 2 

JQ JB 

where the second integral is with respect to the natural surface measure on B. 
It is evident that the forms Q s are monotonic increasing. It may be shown that 
the perturbation term is relatively compact, so that the right-hand side is the 
closed form associated with a certain non- negative self-adjoint operator H s . The 
eigenvalues of H s are increasing real-analytic functions of s by ||. 

Since the quadratic forms are monotonic increasing as a function of s they converge 
to a limit Q defined by 

Q(f) := lim Q s (f) 

s— H-oo 

where we adopt the standard convention [0] that Q(f) = +oo whenever / does not 
lie in the form domain T> of Q. It is clear that 

V ={f 6 V :/+ = /_ on B}. 

This space is exactly W /1,2 (fi), so the operator associated with Q is H := —A 
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acting in L 2 (VL) subject to NBC on dVL. An immediate consequence is that 

lim fi n (s) = n n 
for every n, the limit being monotone. 

Having chosen a suitable increasing sequence = s < si < . . . < s n one ob- 
tains accurate enclosures of the eigenvalues of each H Si+1 from those of H Si by the 
Goerisch-Plum homotopy method. If s n is large enough then the eigenvalues of 
H Sn will be good enough lower bounds of the eigenvalues of H to enable us to 
apply RRTL to obtain accurate enclosures of the eigenvalues of H. Many of the 
details are similar to what we have already discussed, and we concentrate on the 
novelties. 

The upper bounds on all eigenvalues are obtained by RR using a suitable test 
function space lying in the quadratic form domain V or V. An obvious choice 
is to use a finite element subspace in which the elements are linear, but more 
sophisticated elements are probably needed for accurate results. The continuity 
requirement for two elements which have a common edge is suspended if that edge 
lies within B. The restriction of the quadratic form Q s to the test function space 
can be expanded in terms of the values of the elements at the vertices, noting that 
there will be two values at each vertex lying on B, one corresponding to each side 
of B. 

The required lower bounds can only be obtained by TL if one takes a test function 
space lying in the operator domain, which is different for each operator, even 
though every operator H s is equal to —A on its own domain. One can use a 
finite element subspace consisting of C 2 functions, but this has to respect not only 
the Neumann boundary condition on dW but also certain s-dependent internal 
boundary conditions on B. 

Let df± denote the normal derivatives of / on the two sides of B, both taken in 
the direction from the — side of B to the + side. An application of Gauss' theorem 
shows that the internal boundary condition is 

df + (x) = df4x) = s{f + (x)-f4x)} 

for all x G B. 

It is not easy to demonstrate that the above theory works well in practice, without 
a substantial amount of effort writing the relevant code. An appropriate choice of 
the test function space is crucial if one is to get good enclosures, as may be the 
use of a preconditioning procedure, and we leave this to a future publication. 
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